R tiene una inmensa capacidad de graficar y visualizar datos de todo tipo, incluídos datos genéticos.

Las gráficas pueden hacerse desde la base de R (base) o con paquetes especializados en graficar, como lattice, o más recientemente ggplot2 y ggbio. También paquetes especializados en un tipo de datos que incluyen funciones para graficar, como ape para árboles filogenéticos.

En esta sección veremos una introducción a graficar en R usando graphics que es el sistema que viene con base y luego nos enfocaremos en gráficas más complejas y las principales usadas en análisis genéticos. En R se puede hacer mucho más que lo que veremos aquí, recomiendo profundizar.

Una de las mejores formas de aprender a hacer gráficas en R es buscar en internet/libro una gráfica parecida a la que queremos hacer y ver el código. Algunas recomendaciones:

Base Graphics

Estas son las principales funciones para graficar utilizando la base de R. Puedes buscar ayuda de cada una con su nombre, y además en explorar argumentos extras con ?par

Gráficas x,y:**

Dando x, y:

largo<-c(10,20,11,15,16,20)
ancho<-c(5,10,7,8,8,11)
plot(x=largo, y=ancho)

Dando un objeto que tiene dos columnas, se toman automático como x,y:

# ver el contenido de `cars`(una df ejemplo que viene con R)
head(cars)
##   speed dist
## 1     4    2
## 2     4   10
## 3     7    4
## 4     7   22
## 5     8   16
## 6     9   10
plot(cars)

Si queremos especificar qué columnas serán x, y del objeto:

# graficar vel vs distancia
plot(x=cars$speed, y=cars$dist)

Cambiar título de ejes e íconos:

# graficar vel vs distancia
plot(x=cars$speed, y=cars$dist, xlab="Velocidad", ylab="Distancia", cex=0.5, pch=19)

Ejercicio: mira la ayuda de par y explica qué hacen los argumentos cex y pch.

Ejercicio: Repite la figura anterior pero cambiando los puntos por triángulos azules. Necesitarás esto.

Histogramas

Ejemplo con los datos islands (viene con R)

hist(islands)

Barplot

Ahora cargemos un archivo real de datos:

reads<-read.delim("../Practicas/Uni4/data/reads.txt")

Hagamos una gráfica de barras y colorear acorde a info contenida en otra columna:

head(reads)
##   Library    sample  nreads
## 1    Lib1  pobA21_r 1381230
## 2    Lib1    pobB10 1726622
## 3    Lib1    pobB05  819766
## 4    Lib1    pobC05  442508
## 5    Lib1  pobD21_r 1398874
## 6    Lib1 Outa112_r 2024684
barplot(reads$nreads, col=reads$Library)

Definir colores

Los colores que R ocupa para colorear algo están definidos en palette y pueden cambiarse

# Ver colores
palette()
## [1] "black"   "red"     "green3"  "blue"    "cyan"    "magenta" "yellow" 
## [8] "gray"
# Cambiar colores 
palette(c("green", "blue", "red"))

# volver a graficar
barplot(reads$nreads, col=reads$Library)

Además de manualmente, los colores se pueden definir via paletas predeterminadas:

# Cambiar palette a 6 colores del arcoiris
palette(rainbow(6))

# volver a graficar
barplot(reads$nreads, col=reads$Library)

Checa otras palettes parecidas a rainbow en este link, y no te pierdas cómo nombrar muchos otros colores y utilizar otras paletas con más colores en la R Color Reference Sheet

Agregar una leyenda

# Graficar
barplot(reads$nreads, col=reads$Library)
# Agregar leyenda
legend(x="topleft", legend=levels(reads$Library), fill=palette()[1:3])

Nota que legend es una función por si misma (i.e. NO un argumento de plot) que requiere que antes de correrlo se haya corrido plot. Es decir una vez que creamos una gráfica podemos agregar sobre de esta una leyenda. Lo mismo puede hacerse con la función title.

Boxplot

Ejemplo:

boxplot(reads$nreads ~ reads$Library,
        border = c("red", "blue", "darkgreen"))

Graficar árboles filogenéticos

La graficación de árboles filogenéticos se hace con el paquete ape, con el paquete phytools para funcionalidad extendida y con el paquete ggtree. Empezaremos por ape.

Bibliografía recomendada:

Los árboles filogenéticos pueden construirse en R, simularse en R o leerse a R.

Veamos un ejemplo con un árbol simulado:

# Cargar librería
library(ape)
## Warning: package 'ape' was built under R version 3.2.3
# Simular árbol
set.seed(1) # este comando es opcional, sirve para que todas "simulemos los mismos números" y podamos repetir de forma idéntica la simulación cada vez
tree <- rtree(n = 10, rooted=FALSE)

# ¿Qué tipo de objeto es?
class(tree)
## [1] "phylo"
# ¿Qué contiene?
tree
## 
## Phylogenetic tree with 10 tips and 8 internal nodes.
## 
## Tip labels:
##  t10, t6, t9, t1, t2, t7, ...
## 
## Unrooted; includes branch lengths.
str(tree)
## List of 4
##  $ edge       : int [1:17, 1:2] 11 12 13 13 12 11 14 15 15 14 ...
##  $ tip.label  : chr [1:10] "t10" "t6" "t9" "t1" ...
##  $ edge.length: num [1:17] 0.718 0.992 0.38 0.777 0.935 ...
##  $ Nnode      : int 8
##  - attr(*, "class")= chr "phylo"
##  - attr(*, "order")= chr "cladewise"
# Graficar el árbol
plot.phylo(tree, edge.width=2)

La función plot.phylo puede abreviarse como plot. R sabe que debe usar la función plot.phylo y no plot básico (como lo usamos arriba) porque el objeto que le damos es un árbol.

Podemos modificar este árbol de manera similar a cómo lo hicimos en las gráficas anteriores.

Cambiar el tipo de árbol

Ejercicio Revisa la ayuda de plot.phylo y utiliza un argumento de dicha función para graficar el árbol simulado pero que se vea como un abanico y luego como un cladograma (así):

Enraizar el árbol

# plot árbol sin enraizar
plot.phylo(tree, edge.width=2)

# especificar output para enraizar:
tree<-root(tree, outgroup="t2")

# plot árbol enraizado
plot.phylo(tree, edge.width=2)

Chulear el árbol

Ejercicio: Sigue los ejemplos de Phylogenetic trees in R, del blog Sensory Evolution para realizar los siguientes cambios al árbol anterior:

  • Cambia el nombre de las puntas de “t1”, “t2” etc a “especie 1”, “especie 2”, etc y asegúrate de graficar estos nombres en las puntas de tu árbol
  • Agrega un círculo en los nombres de las puntas. El clado de las puntas t9, t6 y t 10 debe tener círculos rosas, la punta t2 gris y el resto verdes.
  • Incremente al grosor de la línea
  • Cambia el color de la línea a verde oscuro

Leer un árbol en R y graficarlo

Podemos leer árboles en formato newick o nexus a R con la función read.tree de ape:

# cargar archivo
maiz.tree<-read.nexus("../Practicas/Uni4/data/tree")

# checar contenido
maiz.tree
## 
## Phylogenetic tree with 165 tips and 163 internal nodes.
## 
## Tip labels:
##  maiz_3, maiz_68, maiz_91, maiz_39, maiz_12, maiz_41, ...
## 
## Unrooted; includes branch lengths.
# graficar
plot(maiz.tree, type="unrooted", edge.width=0.1, cex=0.5)

Vamos a poner colores de acuerdo a la Categoría de Altitud en vez de nombres de muestras.

### Graficar por Categorías Altitud

# leer info extra de las muestras
fullmat<-read.delim("../Practicas/Uni4/meta/maizteocintle_SNP50k_meta_extended.txt")

# ¿Cuántos colores necesito?
numcolsneeded<-length(levels(fullmat$Categ.Altitud))
palette(rainbow(numcolsneeded)) 

# graficar sin nombres de muestras
plot(maiz.tree, type="unrooted", edge.width=0.3, show.tip=FALSE)

# Agregar tip labels que correspondan a las categorías de altitud
tiplabels(pch=20, col=fullmat$Categ.Altitud)
# legend
legend(x= "bottomleft", legend=levels(fullmat$Categ.Altitud), pch=19, col=1:numcolsneeded, cex=1, , bty="n")

Ejercicio: Colorea el árbol anterior por Raza (sin incluir una leyenda porque son demasiadas)

Observa que las muestras 162-165 corresponden a teocintles (Zea m. mexicana y Z. m. parviglumis)

fullmat[162:165, 1:4]
##     OrderColecta NSiembra                     Origen               Raza
## 162           NA       96 Jardin Etnobotanico Oaxaca Zea m. parviglumis
## 163           NA      203              Luis Eguiarte    Zea m. mexicana
## 164           NA      911          Luis Eguiarte 1_1 Zea m. parviglumis
## 165           NA     9107          Luis Eguarte 10_7    Zea m. mexicana

Si quisiéramos colorear todas las tips de maíz, pero no las teocintle podemos especificar esto así:

# graficar sin nombres de muestras
plot(maiz.tree, type="unrooted", edge.width=0.3, show.tip=FALSE)
# Agregar tip labels sólo a los maíces
tiplabels(tip=c(1:161), pch=20, col="black")

Ejercicio Grafica el árbol de maíces de manera que los teocintles sean cuadrados negros y lo s maíces círculos verdes, así:

Mapas en R

En R pueden visualizarse mapas de muchas maneras y de hecho hasta hacer análisis complejos con datos raster, como simulaciones de cambio climático y modelos de distribución potencial de las especies. Aquí sólo cubriremos brevemente cómo graficar un shapefile y agregar puntos.

Carguemos uno de los principales paquetes para manipular mapas:

library(maptools)
## Warning: package 'maptools' was built under R version 3.2.3
## Loading required package: sp
## Warning: package 'sp' was built under R version 3.2.5
## Checking rgeos availability: TRUE

Leer un shapefile

La función readShapePoly de dicho paquete nos permite leer un shapefile de polígonos (para puntos hay que usar otra función parecida ¿cómo crees que se llame?).

Por ejemplo vamos leer a R y graficar el shaphile de las regiones biogeográficas de México:

# leer shapefile
biogeo<-readShapePoly("../Practicas/Uni4/data/rbiog4mgw/rbiog4mgw.shp")

# plot
plot(biogeo)

# colorear por bioregioón
palette("default")
plot(biogeo, border="grey", col=biogeo$PROVINCIA)

## Cambiar colores default a algo más bonito
#¿Cuántos colores necesito?
levels(biogeo$PROVINCIA)
##  [1] "Altiplano Norte (Chihuahuense)"     
##  [2] "Altiplano Sur (Zacatecano-Potosino)"
##  [3] "Baja California"                    
##  [4] "California"                         
##  [5] "Costa del Pacifico"                 
##  [6] "Del Cabo"                           
##  [7] "Depresion del Balsas"               
##  [8] "Eje Volcanico"                      
##  [9] "Golfo de Mexico"                    
## [10] "Los Altos de Chiapas"               
## [11] "Oaxaca"                             
## [12] "Peten"                              
## [13] "Sierra Madre del Sur"               
## [14] "Sierra Madre Occidental"            
## [15] "Sierra Madre Oriental"              
## [16] "Soconusco"                          
## [17] "Sonorense"                          
## [18] "Tamaulipeca"                        
## [19] "Yucatan"
# Generar paleta con colores de RColorBrewer
# ver opciones de colores
library(RColorBrewer)
display.brewer.all()  

# generar paleta
palette(c(brewer.pal(9, "Set1"), brewer.pal(10, "Set3")))

# plot
plot(biogeo, border="grey", col=biogeo$PROVINCIA)
legend("bottomleft", legend=levels(biogeo$PROVINCIA), bty="n", cex=.4, fill=palette())

Ejercicio: cambia el color de las provincias “Tamaulipecas” y “Costa del Pacífico” a otro color.

Agregar puntos a un mapa

Es muy común tener las coordenadas x,y de nuestros puntos de muestreo en un archivo junto con el resto de la info de nuestras muestras. Por ejemplo en el caso de la info de maíces que hemos estado utilizando:

fullmat[1:5,c("Latitud", "Longitud")]
##    Latitud   Longitud
## 1 28.68578 -107.58300
## 2 19.29083  -97.92944
## 3 28.53703 -108.92467
## 4 16.03225  -92.97972
## 5 20.09111 -101.11889

Agregar esta info a un mapa en la misma proyección y sistema de coordenadas puede hacerse con la función points:

# plot map
plot(biogeo, border="grey", col=biogeo$PROVINCIA, lwd=0.8)
# agregar puntos
points(fullmat$Longitud, fullmat$Latitud, pch=19, col="black", cex=0.4)

Ejercicio: Baja un mapa (nivel nacional) que te interese del GeoPortal de la CONABIO, ploteálo y agrega los puntos del muestro de maíz, utilizando una forma de punto diferente para los teocintles, y que los puntos estén coloreados por CategAltitud.

Guardar imágenes:

Método 1:

  • “Abrir un device” con png(), jpg(), pdf() según el formato en que queramos guardar.
  • Hacer la gráfica
  • “Cerrar el device” con dev.off
png("../Practicas/Uni4/out/arbol.png")
plot(maiz.tree, type="unrooted", show.tip=FALSE)
tiplabels(pch=20, col="green3", tip=c(1:161), cex=.5)
tiplabels(pch=15, col="black", tip=c(162:165), cex=.5)
dev.off()
## quartz_off_screen 
##                 2

Método 2:

En R studio darle “Export” en el panel de la imagen y seleccionar un nombre de archivo y demás características.

ggplot2

Las gráficas que hemos visto hasta ahora pueden verse un poco feas de inicio y puede tomar un rato y mucho código arreglarlas a algo hermoso. ggplot2 es un paquete que ahorra este trabajo y que ha comenzado a ser ampliamente adoptado.

Los árboles filogenéticos también pueden graficarse estilo ggplot, con el paquete de Biocounductor ggtree.

No lo cubriremos en este curso pero es bueno que sepan que existe.

Ejercicio

Escribe un script de R para hacer una gráfica parecida a la de la imagen.

Nota que están ordenados con las muestras de maíz de menor a mayor altitud primero y al final (derecha) las muestras de teocintle.

Necesitarás: 1. Asegúrate que tu wd sea Practicas/Uni4/bin, pero NO pongas en tu script setwd.

  1. Los datosQ resultado de correr admixture. Están en Practicas/Uni4/data/maices_admixture.Q". Cada línea corresponde a los valores Q de un individuo de maíz o teocintle en el mismo orden que el archivo "Practicas/Uni4/meta/maizteocintle_SNP50k_meta_extended.txt.
  2. Los datos de Altitud que se encuentran en el archivo maizteocintle_SNP50k_meta_extended.txt.
  3. Para poner en dos páneles uno arriba de otro en una misma gráfica necesitas escribir en una línea de R antes de la línea donde haces tu plot: par(mfrow=c(2,1). Nota: Si quieres volver a tener un sólo panel por gráfica: par(mfrow=c(1,1)